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PRACTICAL PULSE SYNTHESIS VIA THE DISCRETE INVERSE SCATTERING 

TRANSFORM 
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CROSS-REFERENCE TO RELATED APPLICATIONS 

[0002] The present application claims priority to U.S. Provisional Patent Application 

No. 60/433,407, filed December 13, 2002, the contents of which are incorporated herein in 
their entirety. 

FIELD OF THE INVENTION 

[0003] The present invention relates to the generation of pulses using the discrete 

inverse scattering transform (DIST) and, more particularly, to a technique for obtaining RF- 
pulse envelopes that produce a specified magnetization profile, to within a predetermined 
error, for use in MRI and other applications. 

BACKGROUND OF THE INVENTION 
Brief Summary of the Pulse Synthesis Problem 

[0004] The synthesis of RF pulse sequences to produce selective excitation is a 

problem of principal importance in all applications of nuclear magnetic resonance. To discuss 
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this problem, one needs to introduce the Bloch equation without relaxation, which is usually 
written in the form: 



— — = yM x B. 
" (1) 



Here M is the magnetization, B is the applied magnetic field, t is time and y is the 
gyromagnetic ratio. A solution to this equation has constant length which is assumed 
throughout to equal 1. The Bloch equation is usually analyzed in a "rotating reference" frame. 
Ordinarily the rotating reference frame is related to the "laboratory frame" by a time 
dependent orthogonal transformation of the form: 



so that 



cos -sin 0 
sin0(O cos0(*) 0 
0 0 I 



Af(f) = FU)mm. (2) 



One uses m to denote the magnetization in the rotating reference frame. Larmour's Theorem 
implies that if M satisfies (1) then m satisfies 



dm 

- = y m * BcS r (3) 



where 



BefYV) = (f)B(t) + ^«(0 with St(t) = (0, 0 ? 0' (/))'- 



[0005] In most applications of this method, the function 6 is selected to render the z- 

component 

of B e ff independent of time, 



The constant value /is called the offset frequency or resonance offset. If one sets co(f) = 
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co\(t) + ia)2(f) then, in the laboratory frame, the RF-envelope takes the form: 
p{t) - o>U)e?°V\ 

The energy in the RF-envelope is given by 



= j \plO\ 2 d*= j \co(i)\ 2 dt\ 
-*! -«> (6) 

The magnetization profile is a unit 3-vector valued function defined for fe R, 



[0006] In essentially all MR applications m "(/) = [0, 0, 1]' outside of a finite interval. 

In most earlier approaches to pulse synthesis, the flip angle profile, defined as: 

Kfi = cos-Vmf (/)), 

was used to design the pulse, with the phase of transverse magnetization determined, 
implicitly, by the algorithm. Typically, this is referred to as "recovering" the phase. 
[0007] The problem of RF-pulse synthesis is to find a time dependent complex pulse 

envelope co(0 so that, if Z? e ff(/) is given by (5), then the solution of: 



^(/;/) = yOT(/;0^^ff(/:Owith # lim m{f;t) = (0, 0, 1)' 
at t^^x (7j 



satisfies 



lira ET^Owi +/JH2)(S;0,m3(S;0] = l(wf + /«f )C0.i"f C/)l- 

/->«> (g) 

The standard complex notation, mi + imi has been used for the transverse components of the 
magnetization. If ©i(f) + ia>2(f) is supported in the interval [r 0 , *i], then these asymptotic 
conditions are replaced by: 
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The mapping from coi(f) + kd 2 (0 to m°° (as defined in Equations (7) and (8)) is highly 

nonlinear; the problem of pulse synthesis is that of inverting this mapping. 

[0008] To solve the problem of RF-pulse synthesis, it is convenient to introduce the 

spin 

domain formulation of the Bloch equation. Instead of a unit vector m in R 3 , one solves for a 
unit vector ijr'\x\ C 2 . This vector satisfies the 2x2 matrix equation: 

df i 

(10) 

Here go = -[ycdi(0,Ycd 2 (0> v L and o are the Pauli spin matrices: 



a, =[? o]< 01 =[? ol' 03 = [o -l] 

1 J (11) 



Assembling the pieces, it can be seen that ^satisfies: 

?*«»-[-& T]* 8!0 ' (12) 

with 



where q* is the complex conjugate of the complex number q. 

[0009] A simple recipe takes a solution of Equation (12) and produces a solution of 

Equation (7). If ^(£;0 = [^i(£;f)> ^(SiOf satisfies Equation (12), then the 3-vector: 

2 (14) 

satisfies Equation (7). If in addition: 

L J (15) 
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then m satisfies Equation (7). Thus, the RF-pulse synthesis problem is easily translated into 
an inverse scattering problem for Equation (12). This is described in the next section. 
Following the standard practice in inverse scattering, Equation (12) is referred to as the ZS- 
system and q as the potential. 
Scattering theory for the ZS-system 

[0010] Scattering theory for equations like those above relates the behavior of ^(£;0> 

as t — > -~ to that of i/K&\t)\ as t — ► -h». If q has bounded support, then the functions: 

[ e o ]• [A] 

(16) 

are a basis of solutions outside the support of q. If the L'-norm of q is finite, then, it is known 
that there are solutions that are asymptotic to these solutions as t— >±<». 
Theorem 1. If\\q\\i} is finite, then, for every real £ there are unique solutions: 

which satisfy: 

L J L J (18) 



(19) 

The solutions ^i-(£X ife+(€) extend as analytic functions of £ to the upper half plane, {£ : Im£ 
> 0} andjfe-(0> ^i+(£) extend as analytic functions of £ to the lower half plane, {£ : Im£ > 0}. 
The proof of this theorem can be found in Ablowitz, et al., 'The inverse scattering transform- 
Fourier analysis for nonlinear problems," Studies in Applied Math., Vol. 53, (1974), pp. 249- 
315. 

[0011] For real values of £ the solutions normalized at -«> can be expressed in terms of 

the solutions normalized at +«> by linear relations: 

*2-(£; 0 = i+(t; O - «*(*)*2+«; O. (2 0) 
The functions a; b are called the scattering coefficients for the potential The 2x2- 
matrices [^i.(fe-L [^1+^2+] satisfy: 

(21) 

The scattering matrix for the potential <y is defined to be: 
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^-k) -^>J (22) 

It is not difficult to show that: 

It follows from Theorem 1 and Equation (23) that a extends to the upper half plane as an 
analytic function. Also, if {£i,...,£yv} is a list of the zeros of a, then Equation (23) implies 
that for each j there is a nonzero complex number C'j so that: 

= c5* 2+ «» f j = (24) 

The reflection coefficient is defined by: 

a & (25) 

A priori the reflection coefficient is only defined on the real axis. Since |a(£)| 2 + |b(£)| 2 = 1, 
a(£) may be expressed in terms of r(^) as: 

°° log(t + \r(Q\ 2 )dC 



J (26) 

Equation (26) has a well defined limit as £ approaches the real axis. 

[0012] If a has simple zeros at the points {£i,...,£n} (so that a'(^) ^ 0). then the 

norming constants are defined by setting: 



C; = 



J 



(27) 

where the { C )} are defined in Equation (24). The definition needs to be modified if a has 
nonsimple zeros. The reason for replacing {C'j} with {Cj} will become more apparent 
below. The pairs {(£>C/)} are often referred to as the discrete data, or bound state data. 
[0013] The function r(£), for £ e R and the collection of pairs {(£;,C/) : j = 

define the reduced scattering data. Implicitly the reduced scattering data is a function of the 
potential q. In inverse scattering theory, the data {r(0 for £ e K; (£i>Ci )»•••» (£n,C n ) are 
specified, and a potential # is sought that has this reduced scattering data. The map from the 
reduced scattering data, or the mathematically equivalent data, to q is often called the Inverse 
Scattering Transform or 1ST. 

[0014] The RF-pulse synthesis problem will now be rephrased as an inverse 

scattering problem. Since the data for the pulse synthesis problem is the magnetization 
profile m°°, a function of £ = v/2. The solution ifr\. to the ZS-system defines a solution nt\. to: 
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—{v:i) = ym(v:t)xB < :tt{v;0, 

* (28) 

with 

lim »(u;/) = iaO ! l| f , 

i—r —CO 

(29) 

It follows from Equations (19) and (20) that: 

(30) 

Therefore: 

-*■»- terw -™ (31 ) 

[0015] The complex notation for the transverse components of mi. are used. If nt\ 

also satisfies: 

lira le-*""C/»i + f!» 2 )(ip;i), m 3 (»; 01 = Km?* + JjM$°)(ir), mf°(ir)|. 
^co * (32) 

then it follows from Equation (31) and |a(£)| 2 + |b(^)| 2 = 1 that: 



= Mi = iim (^li- -i-/w2i-)(f ^^"'^ 



(33) 



If q has support in the ray (-°°,fi) then: 



is independent of r for f > fi. It is also useful to observe that if r (£) is the reflection 
coefficient, determined by the potential q(t), then e~ 2lx ^r (£) is the reflection coefficient 
determined by the time shifted potential q x (t) =q (t - x). 

[0016] Asw~ (£) is a unit vector valued function, the reflection coefficient r (£) 

uniquely determines m «> (£) and vice-versa. Thus, the RF-pulse synthesis problem can be 
rephrased as the following inverse scattering problem: 

[0017] Find a potential q(t)for the ZS-system so that the reflection coefficient r (Q 

satisfies Equation (33) for all real £ 
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[0018] The pulse synthesis problem makes no reference to the data connected with 

the bound states, i.e. {(£;, C /)}. Indeed these arc free parameters in the RF-pulse synthesis 
problem, making the problem highly underdetermined. 

[0019] The flip angle profile, q>(£); is related to the scattering data by: 

^ =srl (Tw) =2sr,(|6(?,l) - (35) 

SLR pulses are usually designed using the flip angle profile, and the phase of the reflection 
coefficient of the designed pulse is determined indirectly. In this context, the reflection 
coefficient determined by the pulse envelope is an approximation to a function of the form 
e'^^r(Q where r(£) is the "ideal" reflection coefficient. The function (/>(£) is the phase 
error and it is implicitly determined by details of the SLR algorithm. If (pig) is 
approximately linear over the support of r(^) 5 then the magnetization can be approximately 
rephased. In this case, the actual rephasing time comes in part from t u the maximum value of 
t within the effective support of q(t) (i.e., the largest interval where the value of the potential 
is large enough to produce a measurable effect), and in part from the phase factor, e 1 ^^ . 
[0020] A formula for the energy of the pulse envelope may be stated in terms of the 

reduced scattering data. This formula is quite different from the formula taught by Pauly et al. 
that relates the energy in an SLR-pulse to a parameter, ao, which arises as an intermediate 
step in their construction. It does not provide an explicit formula for the energy in terms of 
the magnetization profile and auxiliary parameters. Moreover, it is not evident from this 
formula how the energy depends on the locations of the bound states. The underlying results 
from inverse scattering theory are due to Zakharov, Faddeev and Manakov. 
Theorem 2. Suppose that q(i) is a sufficiently rapidly decaying potential for the ZS-system, 
with reflection coefficient r(0> ond discrete data {(Sj 9 Cj) 9 j = then: 

OO CV ft 

J \q(0\ 2 di = ^ f tog(l +|r(|)| 2 )Jf + 4jrimfe. 

;=l (36) 

Prior approaches to pulse synthesis 

[0021] The oldest method of pulse synthesis is the Fourier transform method. The 

Fourier method is a simple linear approximation for this highly non-linear process. It 
provides usable results for small flip angles, but is largely an ad hoc method. 
Notwithstanding the fact that this method is an approximation for all non-zero flip angles, it 

71 

gives surprisingly good results for flip angles up to about — . Over the past fifteen years two, 

2 
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more exact, approaches to solving this problem have been introduced and explored. The 
slightly older, generally accepted method is the Shinnar-Le Roux- or SL/?-algorithm 
introduced independently by M. Shinnar and his co-workers (USP 5,153,515; USP 5,572,126; 
Shinnar, et al., "The Synthesis of Pulse Sequences Yielding Arbitrary Magnetization 
Vectors," Mag. Res. In Med., Vol. 12 (1989), pp. 74-88; and Shinnar et al. "The Application 
of Spinors to Pulse Synthesis and Analysis," Mag. Res. in Med., Vol. 12 (1989) pp. 93-98) 
and P. Le Roux (USP 5,821,752; P. Le Roux, "Exact Synthesis of Radio Frequency 
Waveforms," Proceedings of 7 th Annual Meeting of SMRM, 1988, p. 1049; J. Pauly et al., 
"Parameter Relations for the Shinnar-Le Roux Selective Excitation Pulse Design Algorithm," 
IEEE Trans. On Med. Imaging, Vol. 10 (1991), pp. 53-65) The second method is to use the 
inverse scattering transform, or 1ST, which was first introduced as a technique for pulse 
synthesis by Grunbaum et al. in "An Exploration of the Invertibility of the Bloch Transform," 
Inverse Problems, Vol. 2 (1986), pp. 75-81, and later refined by J. Carlson, "Exact Solutions 
for Selective-Excitation Pulses," J. of Mag. Res., Vol. 94 (1991), pp. 376-386, and "Exact 
Solutions for Selective-Excitation Pulses. II. Excitation Pulses with Phase Control," J. of 
Mag. Res., Vol. 97 (1992), pp. 65-78; and Rourke et al., 'The Inverse Scattering Transform 
and its Use in the Exact Inversion of the Bloch Equation for Noninteracting Spins," J. of 
Mag. Res., Vol. 99 (1992), pp. 118-138. There are also approaches which use optimization 
techniques from control theory, but such approaches are not particularly relevant to the 
subject matter of the invention and thus will not be discussed further. 

[0022] In the SLR-method, one approximates not the desired magnetization profile, 

but 

the flip angle profile derived from it. It is approximated by a rational function, r 0 of e if ^^ . 
Using standard polynomial design techniques, one finds a hard pulse, that is, a sequence of 
equally spaced 5-pulses: 

at 

which produces an approximation, mo, to the desired magnetization. Note, however, that only 
the flip angle is directly controlled. The phase of the excitation is determined implicitly by 
the algorithm used to construct r 0 . By specifying only the flip angle profile, the SLR 
approach retains more direct control over the total duration of the RF-pulse. 
[0023] The forward scattering analysis has a well-defined limit for potentials of the 

form q 0 (t). In this case, the scattering matrix takes the special form: 
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where A 0 (z) and S 0 (z) are polynomials of degree N - I. The reflection coefficient r 0 (0. 
defined by </o(0> is the periodic function of period A' l 2n given by: 

e-™*S Bale- 1 **) 



ro«) = 



4o(«-'*) (3g) 



[0024] Thus, the SLR algorithm has two parts: 1. Find polynomials of the given 

degrees, so that ro(0 is, in a certain sense, an approximation to 2. Use a recursive 
method for determining the coefficients, so that qo(t) given as: 

A T 

i-i (39) 

has reflection coefficient ro(0. 

[0025] Of course, a sum of 5-pulses is nonphysical, requiring infinite energy to 

realize. The RF-envelope that is actually used is a "softened" version of qo(t). For example, 
one could replace each \ij6(t - jA) by a boxcar pulse of width A with the same area, leading to 
the softened pulse: 

A' 

im| (40) 
While the difference #o(0 - <?i(t) can only be made small in the sense of generalized 
functions, the difference | ro(£) - | can be made pointwise small, for £ in a fixed interval, 
provided that none of the is too large. This is what is usually meant by the "hard pulse 
approximation" described by Shinnar et al. 

[0026] Careful control of the duration is important in applications that involve 

imaging samples with a short F2, or spin-spin relaxation time. Spin-spin relaxation only 
becomes important once the magnetization has a significant transverse component. At least 
for minimum energy pulses, this occurs near the peak of the pulse. The duration of the pulse 
after the peak is almost equal to the rephasing time. Hence, in order to design pulses for 
experiments with a short T2, it is really only necessary to control the rephasing time, which is 
possible in 1ST pulse design. In practice, a pulse designed using the 1ST method has only 
slightly longer duration than the SLR pulse, with the same design parameters. In numerical 
simulations, an 1ST pulse appears to be slightly less susceptible to errors caused by spin-spin 
relaxation, than the comparable SLR pulse. 
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[0027] The SLR-method does not directly control the phase of the transverse 

magnetization profile nor do the "auxiliary parameters," inherent in the pulse design problem, 
arise in this approach. In truth, the SLR approach does not solve the pulse design problem, 
per se, though in many important instances, it produces a result which is adequate for most 
current applications. 

[0028] From a practical perspective, both 1ST and SLR require considerable 

computation. From the inverse scattering formalism, it is clear that there are many possible 
RF-envelopes that will produce the same magnetization profile. There are also auxiliary 
parameters in the SLR method. Changing the auxiliary parameters that appear in SLR 
changes the phase of the magnetization profile, only leaving the flip angle profile unchanged. 
So changing the auxiliary parameters in SLR actually produces solutions to different pulse 
design problems, whereas varying the auxiliary parameters in 1ST produces different 
solutions to the same pulse design problem. In either case, some additional criteria are needed 
to determine how the auxiliary parameters should be chosen. The "real" problems of pulse 
synthesis are therefore: 

(a) To choose the auxiliary parameters leading to an RF-envelope, "optimal" for a 
given application, which produces an approximation to a given magnetization profile; and 

(b) To provide an algorithm that produces the optimal pulse. 

Point (a) entails having an algorithm that can find the pulse, with a given magnetization 
profile, and an arbitrary specification of the auxiliary parameters. Prior to the work of the 
present inventors, this problem was unsolved in the MR literature. This "real" problem of 
pulse synthesis is solved herein and used to synthesize pulses of the type usable for MRI, for 
example. 

[0029] In view of the above, there is a need in the art for a pulse generation algorithm 

that is as efficient as SLR though more flexible and that provides more direct control over the 
final result in terms of the input data. 

SUMMARY OF THE INVENTION 

[0030] The above-mentioned needs in the pulse synthesis art are addressed by the 

techniques of the invention. The techniques described herein show the following: 
(a) how to construct pulses with a specified magnetization profile and rephasing time. 
This is in contrast to SLR where only the flip angle profile is specified, and the phase of the 
transverse magnetization is "recovered." 
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(b) how to encode the auxiliary parameters so as to regain the full space of solutions to 
the pulse design problem, known to exist by the inverse scattering approach. The problem of 
pulse design is thus reduced to solving a family of integral equations known in the 
mathematics literature as the Marchenko equations, 

(c) for a given magnetization profile, how to choose the "auxiliary" parameters in the 1ST 
approach to obtain the RF-envelope with the minimum possible energy requirement. 
[0031] The left Marchenko equation is also introduced to address the numerical 
difficulties caused by usage of auxiliary parameters. Two different, but equivalent 
formulations of the inverse scattering problem, are used. This leads to the left and right 
Marchenko equations. Previously, only the right Marchenko equation had been applied to the 
problem of RF-pulse design. 

[0032] Also, the discrete inverse scattering transform (DIST) is explained herein, 

which provides the unique solution to a discretization of the Marchenko equations, with 
arbitrarily specified magnetization profile, and auxiliary parameters. The algorithm is quite 
stable, and as efficient as the standard (very efficient) approach used in SLR. The DIST is a 
particular discrete approximation of the full inverse scattering transform. 
[0033] Two approaches to solving the pulse design problem are proposed in 

accordance with the invention. The first approach uses a straightforward discretization of the 
continuum model. The second approach uses the hard pulse approximation, as is also done in 
SLR. 

[0034] The present invention differs from the SLR approach in that the algorithm of 

the invention approximates the magnetization profile directly, and not the flip angle profile. 
The first approach uses a finite sum approximation to the Marchenko equations, while the 
second approach uses a recursive algorithm that is a different discretization of the Marchenko 
equation formalism. Both algorithms use the left Marchenko equation to deal with significant 
numerical instabilities that arise in the presence of bounded states if only the right equation is 
used. Using the Marchenko equation formalism is believed by the present inventors to be 
important if one wishes to freely specify the auxiliary parameters, and thereby obtain the full 
set of solutions to the pulse synthesis problem. 

[0035] The methods described herein also provide a way to retrospectively analyze 

existing pulses to see if their energy consumption can be reduced. A variety of tools are 
provided for adapting the RF-envelope for specific purposes. For example, the method of the 
invention gives a method to reduce the rephasing time of the pulse, at the expense of a known 
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increase in energy. The method of the invention also provides a method to study the 
dependence of the magnetization profile on the amplitude of the RF-envelope. 
[0036] In a presently preferred application of the invention, RF pulses are generated 

using the techniques of the invention for magnetic resonance imaging and spectroscopy. 
However, as more sophisticated applications of magnetic resonance are found (e.g. hybrid 
forms of spatial encoding), and attempts are made to use less homogeneous background 
fields, the more robust and direct approach to pulse synthesis in accordance with the 
invention will be preferred to the existing methods. Indeed, these more demanding 
applications may well require precise control of the phase of the magnetization profile. Such 
control is attainable using the method of the invention. The analysis of, and solutions to the 
pulse design, provided by the present invention also apply to any two-level quantum system 
or subsystem described by the spin domain Bloch equation, such as systems in the emerging 
field of spintronics and such potential applications as quantum computing where it is also 
very important to control the phase of the excitation. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0037] The present invention will be apparent from the following detailed description 

of the invention taken in conjunction with the accompanying drawings, of which: 
[0038] Figures l(a)-6(f) illustrate various truncations of a 90% minimum energy, 1ST 

pulse and the magnetization profiles they produced. 

[0039] Figure 2 shows a comparison of in-slice versus out-of-slice errors for two 

methods of pulse design for two different choices of transition width and rephasing time. 
[0040] Figure 3(a) illustrates a 90° pulse with 2ms rephasing time, 8 2 = 0.1 and the 

transition width is 0.2 KHz. 

[0041] Figure 3(b) illustrates SLR versus 1ST for a 90° pulse with 4ms rephasing 

time, 82 = 0.01 and the transition width is 0.15 KHz. 

[0042] Figure 4(a) illustrates T 2 = 10 ms for a 90° pulse with 2ms rephasing time, 62 = 

0.1 and the transition width is 0.2 KHz. 

[0043] Figure 4(b) illustrates T 2 = 10 ms for a 90° pulse with 4ms rephasing time, 62 = 

0.01 and the transition width is 0.15 KHz. 

[0044] Figure 5(a) illustrates a minimum energy, selective refocusing pulse designed to 

flip the equilibrium magnetization through 179° using the techniques of the invention. 
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[0045] Figures 5(b)-5(d) illustrate the effect on the pulse in Figure 5(a) on the initial 

vector [1,0,0], [0,1,0], and [0,0,1], respectively. The x -components are shown with a solid line, 

the y-components with a dashed line, and the z-components with a dotted line. 

[0046] Figure 6 illustrates a variety of self refocused 90° pulses that were designed by 

using different rational approximations, r, for the desired reflection coefficient. 

[0047] Figure 7(a) illustrates the minimum energy pulse with the given magnetization 

profile, M. 

[0048] Figure 7(b) illustrates the pulse with magnetization profile, M, and a bound 

state at 0.5/, with norming constant 1.0. 

[0049] Figure 7(c) illustrates the pulse with magnetization profile, M, and bound states 

at 0.5/ + 1.4, and L0i~L0, with norming constants 1, and -10. 

[0050] Figure 7(d) illustrates the pulse with magnetization profile, M, and bound states 

and i-1, i, and i+1, with norming constants 1, 2 and 1. 

[0051] Figure 8(a) shows the minimum energy 1ST pulse for the Mexican hat 

example. 

[0052] Figure 8(b) shows the transverse magnetization defined by jpfl;) for the 

minimum energy 1ST pulse of Figure 8(a). 

DETAILED DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS 
[0053] A detailed description of illustrative embodiments of the present invention will 

now be described with reference to Figures 1-8. Although this description provides detailed 
examples of possible implementations of the present invention, it should be noted that these 
details are intended to be exemplary and in no way delimit the scope of the invention. 
[0054] It also should be understood that the methodology disclosed below can be 

implemented with the aid of software development tools that execute on a computing device. 
Such software development tools may be embodied in the form of program code (i.e., 
instructions) stored on a computer-readable medium, such as a magnetic, electrical, or optical 
storage medium, including without limitation a floppy diskette, CD-ROM, CD-RW, DVD- 
ROM, DVD-RAM, magnetic tape, flash memory, hard disk drive, or any other machine- 
readable storage medium, wherein, when the program code is loaded into and executed by a 
machine, such as a computer, the machine becomes an apparatus for practicing the invention. 
The present invention may also be embodied in the form of program code that is transmitted 
over some transmission medium, such as over electrical wiring or cabling, through fiber 
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optics, over a network, including the Internet or an intranet, or via any other form of 
transmission, wherein, when the program code is received and loaded into and executed by a 
machine, such as a computer, the machine becomes an apparatus for practicing the invention. 
When implemented on a general-purpose processor, the program code combines with the 
processor to provide a unique apparatus that operates analogously to specific logic circuits. 
OVERVIEW 

[0055] The present invention relates to a practical implementation of a method for 

selective RF-pulse synthesis by solving the inverse scattering problem. 
[0056] In particular, given the reduced scattering data {r (£), 1; (£y, C7),. . ., (<f 

Cn)}, with <f , £ £ j , the following function may be defined: 

| ? A' 

" - 03 /ml (41) 

This is the inverse Fourier transform of r (£) with a correction added to account for the bound 
states. The finite sum in Equation (41) is exponentially decreasing as t tends to + <». For 
each fe R,; the operator F t may be defined by: 

oc 

Fth (s) = f f(s + y)h{y)dy for s € [i , oo). 

(42) 

One may denote the adjoint of F t , as an operator on L 2 ([t, «.)), by F*. For each t e 
(_oo 5 oo) the (right) Marchenko equation is the integral equation, for a function k , (s) 
defined for s e [r, ~), given by: 

Kid +f;*F l )* l k*> = r u + *>> ^ 43) 

or more explicitly: 



cor co ~~i 



(44) 

The solution of the inverse scattering problem is given in the following theorem: 
Theorem 3. Given a smooth, rapidly decaying reflection coefficient r (£) 5 and a finite set of 
pairs {(£ j 9 C j) :j = 1,..., N], with the {£ j) distinct, Im £j > 0 and Cj£0,for j = 1,..., Af, 
Equation (43) ftas a unique solution for every re K. //• 

</(/) = -2*, (/), (45) 

tfien rfte ZS-system, with this potential, has reflection coefficient r. It has exactly N bound 
states at frequencies . ., £v}, and tfze relations of Equation (24) Zio/d ar //tere points. 
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[0057] As noted above, r = b/a\ given r; and the locations of the bound states, 

Equation (26). May be used to determine a, and therefore b. For ^e8, define: 

««> (46) 

and let . be the zeros of a in the upper half plane. The kernel function for the left 
Marchenko equation is defined by: 



CO 



5(0 = ^: j r^e^dt-if^Cje^^ 

(47) 



where the coefficients * c '' are given by: 

t 



(48) 

The finite sum of exponentials in Equation (47) decays as t tends to -<». The values {a* (5;)} 
are easily computed using Equation (26), implying that: 



„. , _J i-r ( tk-j} \ F i 7 to8(H-|rg)| 2 )dc l 



(49) 

If the bound states are close together, then | a v (€*) 1 2 may become quite small; this, in turn, 

makes c j quite large. So even though the exponentials in Equation (47) decay rapidly as t — ► 

-°o the finite sum in Equation (47) can be very large for / near to zero. 

[0058] For each t 6 K define an operator from L 2 ((-<~,f]) to itself by: 

> 

G,I(s)= j g(s + y)Hy)dy. 

-~ (50) 

For each t e M , the left Marchenko equation is: 

l(ld+C,G;)/,l<s)=-g(s+0, fbrse (-00,/J, 

or more explicitly: 

-* 0, L-oc. J (52) 

If l t (s) solves this equation, then the potential defined by: 

(53) 

has reflection coefficient r(£)> It has exactly N bound states located at frequencies ..,£yvh 
and the relations in Equation (24) hold at these points. 
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[0059] The input to the 1ST pulse design process in accordance with the invention is 

an ideal magnetization profile, m~(v) which, via Equation (33), defines an ideal reflection 
coefficient r,(Q. Frequently, the ideal reflection coefficient is a discontinuous function. The 
following five step procedure will provide a general class of solutions to the inverse 
scattering problem. 

1. The first step in 1ST pulse synthesis is to choose a smooth function with bounded 
support, r d (g) which is an approximation to r,-(5). So far as the algorithm is concerned, the 
choice of approximation, r<$0, is constrained only by the need to have its Fourier transform, 

decay reasonably quickly. There is no need to use any special type of function such as a 
rational or meromorphic function. As with SLR pulse design the "transition width" and ripple 
directly effect the duration of the final pulse. 

2. The next step is to choose the locations of the bound states and the norming constants, 
{(£i,C,»i = 1,...,7V}. For a minimum energy pulse, no bound states are used. For a self 
refocused pulse, r d (Q is usually taken to be a rational function. The are then the poles of 
r d (g) in the upper half plane, and the {C, } are the residues of r rf (£) at these poles. 

3. Given the data { r rf (0; C\ ),..., (£yv,Cyv)}, the left reflection coefficients and left 
norming constants are determined, using Equations (46) and (48), respectively. With this 
data, one can determine the kernel functions, /and g; for the right and left Marchenko 
equations, using Equations (41) and (47), respectively. In practice, these functions are 
computed at sample points by using the Finite Fourier Transform (FFT). The functions Mf(i) 
and M g (t) are computed using, for example the functions: 



Af/(/)=sup|/(s)| t //(/) = f \f(s)\ds 

t 

t 

Af v (0=sup|g(A% /„(/) = f \g(s)\ds. 



(54) 

and the "switching" time t 0 is determined as the unique solution to the equation M/(t 0 )= 
M,(t 0 ). 

4. The left and right Marchenko equations are approximately solved to determine 
{q d (tj)} at sample points (tj) belonging to a sufficiently large interval [t\t + ]. This interval 
should contain the effective support of q d (t) which can be determined, a priori, by using the 
estimates: 



2Af/(2/) 2iW,,(2#) 



(55) 
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The right equation is used for tj> x 0 and the left equation for r 7 < x 0 . The fineness of the 
sample spacing, A, is determined by the range of offset frequencies present in the sample. As 
is well known, using a sum of boxcar pulses of duration A to approximate qjjt) produces 
sidelobes in the magnetization profile at intervals of length 2A" 1 7c. 

5. The duration of the pulse produced by the 1ST (i.e., the effective support of q^f)) may 
exceed what is allowable in the given application. In this case the pulse is truncated, and 
usually windowed with a smoothing function, e.g. cos 2 (B(t-t 0 )). This leads to changes in the 
excitation profile, usually a decrease in selectivity as well as phase errors, as shown in Figure 
1. Figures l(a)-l(c) illustrate various truncations of a 90% minimum energy, 1ST pulse and 
the magnetization profiles produced for sharp truncation by .001% criterion (Figure 1(a)), 
sharp truncation by 1% criterion (Figure 1(b)), and smooth truncation by l%n criterion 
(Figure 1(c)). The transverse magnetizations produced by the truncations of Figures l(a)-(c)) 
are illustrated in Figures l(d)-l(f), respectively. Sometimes a better result is obtained by 
starting over, replacing r d (Q by a "smoother" function, i.e. a function with a wider transition 
region. 

[0060] The approach of the invention is a particular discrete approximation of the full 

inverse scattering transform (1ST) approach described above which is, in a sense, a full 
implementation of the hard pulse approximation within the inverse scattering approach. In the 
discrete inverse scattering transform (DIST) approach of the invention, the choice of r<*(£) is 
further constrained to be of the form: 

(56) 

As in the SLR approach, the present invention uses a hard pulse approximation, directly 
designing an RF-envelope of the form: 

CO 

(57) 

Each algorithm is then a two step process: 

Step 1: Approximate the "ideal" magnetization profile by a periodic function from a 
specified class. 

Step 2: Using the result of step 1, find a hard pulse which produces the given 
approximation to the magnetization profile. 

[0061] In practical applications, a "softened" version of qci is actually implemented. For 

example, one might actually use: 
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00 

In order for the softened pulse to produce a magnetization profile close to that produced by qo, it 
is necessary that the coefficients {cojj be uniformly small. As these coefficients are O(A), this 
reduces to the requirement that the time step, A, be sufficiently small. This is what is often 
called the hard pulse approximation. 

[0062] The scattering formalism used to analyze the spin domain Bloch equation, with 

a reasonably smooth potential, has a well defined limit for a potential of the form given in 
Equation (57). One obtains a scattering operator, which can be written in the form: 

™ W -L*« •« J' (58) 

The physical frequency, z, is twice the spin domain frequency, i.e., z=2£. For a hard pulse, 
the coefficients, a(z) and b(z) are periodic functions of period A'^tt, and therefore can be 
expressed as functions of w = c lzA . 

[0063] The first step in the SLR approach is to find a polynomial in w that 

approximates the flip angle profile, essentially \b(w)\. A common technique is to use the 
Chebyshev or equiripple approximant, found using a Remez algorithm. This then leads to 
polynomials (a(w), b(w)). In accordance with the invention, the reflection coefficient defined 
to be: 

_ b 

(59) 

and is used directly. The magnetization profile M(z) = [Mx(z), M y (z), M z (z)J' determines, and 
is determined by r as follows: 

w 



LT7FFJ 



M{z) = 

(60) 

[0064] For step 1, a Remez type algorithm may be used to find a Chebyshev 

approximation for r. This is a very important difference between the approach of the invention 
and the SLR approach. By approximating r, rather than \b\, the technique of the invention 
retains control on the phase of the final magnetization profile. To some extent, the technique 
of the invention sacrifices direct control on the duration of the pulse. To be more precise: 
while the full duration of the pulse is not specified in advance, the technique of the invention 
does control the rephasing time of the designed pulse. In practice, the rephasing time is usually 
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close to the length of the part of the pulse after the peak. It is during this part of the pulse that 
there is a significant transverse component, which is, in turn, is subject to the spin-spin (T 2 ) 
relaxation process. In practice, the pulses designed using the 1ST approach have only a slightly 
longer duration than the SLR pulse with comparable design parameters. Generally, the 1ST 
pulse produces a markedly cleaner, and more accurate magnetization profile. 
[0065] In the SLR approach, a very efficient recursive algorithm is used to determine 

the hard pulse from the polynomial approximations to a and b. However, there are, in fact, 
many choices for (a, b) y producing different pulses, which nonetheless generate the same 
magnetization profile. In the inverse scattering language the additional data are called "bound 
states." The bound states arise only implicitly in the SLR formalism. Moreover, in the SLR 
approach, the introduction of bound states changes the phase of the transverse magnetization, 
preserving only the flip angle profile. 

[0066] Accordingly, as step 2 of the algorithm of the invention, a different, very 

efficient recursion is used that generates the pulse from the reflection coefficient and the data 
connected to the bound states. In this approach, the bound states are auxiliary data, which do 
not change the magnetization profile. They capture the true nonuniqueness of the solution to the 
pulse design problem. The recursion of the invention is philosophically quite different from the 
SLR recursion, and is, instead, inspired by the Marchenko equation formalism used in 
continuum inverse scattering theory. 

[0067] The introduction of bound states often leads to ill conditioned problems, with 

exponentially growing condition number. An important aspect of the approach of the invention 
is that it uses both the left and right Marchenko equations. While ill conditioned systems still 
arise, they remain amenable to numerical solution. 

[0068] To summarize, the approach of the invention differs from SLR in three basic 

ways: 

(a) The invention obtains an approximation to the magnetization profile, rather than the flip 
angle profile. 

(b) The invention does not directly control the duration of the pulse but does control the 
rephasing time, and this is adequate for most applications. In practice, the 1ST pulses 
have durations that are comparable to SLR pulses with similar design parameters. 

(c) The invention uses a different recursion to construct the pulse from the scattering data, 
that is r and the bound states. This recursion allows for the arbitrary specification of 
the data connected to bound states. 
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[0069] The present patent application is organized as follows. The next section 

introduces concepts and notation connected to the analysis of hard pulses. The following section 
explains what sort of magnetization profiles can be obtained using hard pulses under various 
constraints. The approach to approximating the magnetization profile is then explained. In the 
following sections, the derivation of the discrete version of the left and right Marchenko 
equations is explained, and the recursive algorithms used to construct the pulse are explained. 
After deriving the algorithm, a variety of examples are given. 
PRELIMINARY DEFINITIONS AND NOTATIONS 

[0070] This section presents some notation used herein. Though the Bloch equation is 

not discussed explicitly, all calculations herein are done implicitly in a rotating reference frame. 
[0071] Precession and rotation operators. If one considers a magnetization state 94: 

M - if to be a unit vector valued function of the offset frequency z e R (measured in radians 
per unit of time) with values in i? 3 , then 94 denotes the vector space of all such functions. SH 
denotes the space of smooth functions from R to SO (3), the group of rigid rotations. 
[0072] An element O(z) of 91 acts on an element M(z) of 94 by the rule: 

O -M(z) = 0(z)M(z). 

[0073] Two families of elements of 9i may be defined For t e R, let P t : fW- 94 denote 

free precession for t units of time: 



St? 



= | SN^"' 2 0 

0 1 



(61) 
For coe Q let 

denote a rotation of \a\ radians around the [SRR.'SR.* 0] -axis. For simplicity, a matrix representation 
for is not given at this stage, but unlike P h the operator Rq> does not have a z-dependence. 
[0074] Scattering theory. A hard pulse Q : 2 - C is defined to be a doubly infinite 

sequence of complex numbers. If one is interested in a pulse of finite duration, then it is required 
that coj = Q(j) be non-zero for only finitely many values of J. In this case, the index of the first 
entry after the last non-zero entry signifies the number of rephasing time steps which should be 
performed immediately following the application of the pulse. More explicitly, suppose that o>p_ 
N ,... s C0p.i are the non-zero entries of Q. Then Q is a hard pulse of duration N with p rephasing time 
steps. After rephasing, the equilibrium magnetization state is rotated to 
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(62) 

Here A > 0 is the time step, and P. pA is the rephasing operation. The function Mq is called the 
magnetization profile for Q. 

[0075] For theoretical purposes, it is useful to define Afn for certain hard pulses, which 

are not necessarily of finite duration. A hard pulse Q is admissible if, for each j e there exist 
elements of SH, r+j(co), r.j(co), with certain analyticity properties that satisfy the recursion: 

T±j+i = PaRvjT+j , (63) 

and normalized by : 

lim r+;/l_; A = Id. 

j~±oc J (64 ) 
Finite duration pulses are always admissible. Generally speaking a pulse is admissible if it decays 
rapidly enough as j- ±°°. 
[0076] For admissible £1: 

Sq r= r~ jT-j f»r any JgZ, 

and 



[!] 



The express ton T + j T-j is independent of j because 

If coj = Q(j) vanishes outside the interval [p — TV, p), as above, then Equation (64) implies 
that T+j = P jA for all j > p and that T.j = P, A for all j <p — N. Therefore, Sq = P. p aT. >p . 
Now using the recursion (15): 

Sq = P-p&P&Rci>p-\ ' ' 'PaRq>p-N+\P&Rq>p-nP(p-N)A> 

which implies that the new definition of Mq is consistent with the original definition in the 
finite duration case (note that P( P -n>a fixes the equilibrium magnetization). 
[0077] The operator Sq is called the scattering operator, and the map CI h-> Sq is called 

the scattering transformation. Since the scattering operator specifies the magnetization profile 
as a function of the offset frequency, the problem of selective excitation pulse synthesis is the 
problem of inverting the scattering transformation. There is a direct analogy between the 
matrix valued functions [T±j( a ))^ and the normalized solutions, y/ ±(z; t), to the spin domain 
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Bloch equations. The scattering operator defined above is then a direct analogue of the 
scattering operator for the Bloch equation, which is defined as: 

s(z) = y/ + (z;ty ] t). 

This operator is also independent of t. 

[0078] The projection operators n + and EL. Define w :M - S l by w(z) = e izA y so 

that every function fin L 2 (S J ) can be written in terms of its Fourier series: 

CO 

The notation J?~(f)(n) may be used to denote /„, the nth Fourier coefficient of /. 

[0079] Let n + : L 2 (S l ) - ht(S l ) and n_ : L 2 (S l ) - H~(S l ) be the projections defined 

by: 



n + ( £ ■»"> = X>«* n 

a=— CO r?—\ 
CO -t 

n-( £ f a w") = Y, f*^' 



and let fi + : L *< 5, > - " +(5 '> - d ft " : ^ "* * ~ (5,) be the projections defined by 



CO | CO 



2 

CO i —J 

ft-c 5Z ■ / '» u ' c > = 2 /o+ ^ 

Notice that H . c .(\*£ } consists of the functions which extend to be analytic function of (O in 
the unit »^D. 

[0080] The operator n c . Suppose that F is meromorphic in U cr C, with a pole at £ 0 - 

The function has a unique representation of the form: 

A* 



(65) 



Here ft( <2>J is a meromorphic function in U, which is analytic in a neighborhood of Co. > 
denotes the finite sum of singular terms in Equation (65). Suppose now that f(co) 
is a function in L 2 (S X ) with a meromorphic extension F to the unit disk. If £ = 
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{£i,...,£m} is a subset of • then the function F has a unique representation of the 
form: 

m 

7=1 (66) 

with h(co) is a meromorphic function in ^ which is analytic in a neighborhood of each Q. n c 
/e L 2 (S') is defined by: 



(67) 

If C is written as the disjoint union of two subsets C then: 

n c f=u c f+u c f. 

Using the Laurent expansion it is easy to see that Tl^f = njl^/. 
[0081] The following lemma is used several times: 

Lemma 1. Let /e L 2 (S J ) have a meromorphic extension, F, to the unit disk, and let C,be a 

finite subset of B * If the poles ofF lie in the set £ then U c f= U_f Ifg e fr+ (S J ), 
then 

n c (f g ) = n_[g(n c f)i (68) 

Proof The first claim follows easily from Equation (67), and the Laurent expansion for (w 
— for j e N. For the second claim, it suffices, by linearity, to consider F(w) = (co — 

Let G denote the holomorphic extension of g to °. At £o one gets the representation: 
where H is holomorphic in the disk. Then: 

This clearly equals Tl_[g(Xl^f)]. 

REALIZABLE MAGNETIZATION PROFILES 

[0082] This section considers which magnetization profiles are realizable with a hard 

pulse. This problem has been solved by Shinnar et al. in "Inversion of the Bloch Equation," J. 
Chem. Phys., Vol. 98 (1993), pp. 6121-6128, for hard pulses of fixed finite duration. 
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[0083] First, the Shinnar-Leigh result describing the magnetizations that are attainable 

using a finite duration pulse are provided. Next, pulses of infinite duration, but finite rephasing 
time are considered, and the reflection coefficient is introduced. Finally, infinite duration pulses 
with infinite rephasing time are discussed. 

[0084] Finite duration pulses. Parts of this theorem are proved below within the 

framework of inverse scattering. 

Theorem 4, Let Q be a finite duration pulse. Then there exist functions: 

] ' l ' p (69, 70) 

a 0 >0 (71) 
\a\ 2 + \b\ 2 = 1 on the unit iv-circle t (72) 



with 



and 



such that 



Ssi = 



-3Krf l ) 2 + lr 2 .| ma'f-tfl 23\a'b j 

-ZSab -2Wab |a| 2 - 1&| 2 J 

?ma*b 

\a\ 2 -\b\ 2 J 



Conversely, if a and b are any functions of the form of Equations (69) and (70), satisfying 
Equations (71) and (72), then there exists a pulse £1 of duration N, and rephasing time p such 
that Sn and Mo are given by Equations (73) and (74). 

[0085] In the spin domain formulation, the action of the real 3x3 matrix in Equation (73) 



corresponds to the conjugation 
[0086] Infinite duration pulses with finite rephasing. Let Q : 2 ~* c be a hard pulse, and 

let p G Z ' As before, Q. has p rephasing time steps if G) y = f2(j) vanishes for all 7 > p. The following 
theorem provides an answer to the question of which profiles Mq are realizable by a pulse with p 
rephasing steps, in terms of the data (a, b). A more useful answer is given in Theorem 6. 
Theorem 5. Suppose that 
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a - Yl 0 ^' 

;=o 



(75) 



CO 



(76) 

2/ C K l„|2 . Iul2 



are functions in L (S ) satisfying |a| + |b| =1 and ao > 0. It is also assumed that a has finitely 

i 

many zeros in the unit disk, that a and b have no common zeros in the unit disk ,and that ~ is 
bounded on S 1 . Then, there exists a hard pulse Q with p rephasing time steps such that 



S n = 



-3| (a* ) 2 + & 2 | S»I(fl*) 2 - & 2 J 2*K«*6 
-2£ta& -23ta& \a\ 2 -\b\ 2 

22a* b 1 
|n| 2 - |fc| 2 J 



[0087] 



The Reflection Coefficient Let a and b be as in Theorem 5, and let Q be the 



corresponding hard pulse. The function r = a . is called the reflection coefficient. On the 

|a| 2 = ' and|&| 2 =-±k. 
unit circle, From this one can compute: 



2a* b = 2|«| 2 - = 
a 

11 11 T+FF 



2r 



7- 



Therefore Mn can be written in terms of r as: 



1-p 



[0088] In fact one can specify the single function r, for a and can be recovered from the 

reflection coefficient, and the locations and multiplicities of the zeros of a. Let £l,...,£m e B 
be the zeros of a in the unit disk with multiplicities <ij,...,<i m e N, and set: 



so that a is analytic and non-vanishing in D - This implies that log a. is analytic in °* On the 
unit circle: 

|«|=|<7| = (l+|r|-)- ,/2 

taglfil = -ylog(l+|r| 2 ). 
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Because a is analytic: 



logo = -n + lpg(l +|r|-). 
Therefore, a can be expressed in terms of r: 



Then ft can be obtained as the product: 

b = ar, (78) 
[0089] If a vanishes at the points {£i,...£, n } then, by analogy with the continuum 

case, there are bound states located at these points. 

[0090] If one starts with an arbitrary reflection coefficient re T (S 1 ) and a finite 

set of points {Q} c D * and positive integers fdjj, then there is no reason to expect the 
resulting function b to be of the form of Equation (76). The following Theorem gives 
criteria, which the reflection coefficient should satisfy, in order to require a given number 
of rephasing time steps. 

Theorem 6. Consider the data (r; £i,...,£ m ; di,...,d m ) forre L°°(S l ) } ^e ° \ {0}, and 
dj e N. Let: 

;=o 

b = f] bpJ 

be the functions obtained using Equations (77) and (78). Then a and b satisfy the hypothesis of 
Theorem 5 if and only if rziP~ l is meromorphic in the unit disk with poles at Ci>--->Cm and 
multiplicities di,...,d m . 

Proof. Suppose that a and b satisfy the hypotheses of Theorem 5. Then bj = 0 for all j < 

-p, so that bw^' is analytic in the unit disk. Therefore, . . . . a is meromorphic in 

. 1 
the unit disk. Clearly, rvtA has the same poles and multiplicities as a since a and b are 

assumed to have no common zeros. 

[0091] Conversely suppose that rvf' 1 is meromorphic in the unit disk with the 

same poles and multiplicities as a . Then bwf' = arvf~ is clearly analytic in the unit disk, 
which implies that bj = 0 for all j < — p. The other hypotheses of Theorem 5 are easily 
verified. 
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[0092] Theorem 5 may now be expressed in terms of the reflection coefficient: 

Suppose that r is a bounded function in L^S 1 ), and suppose that xziP A is meromorphic in the 
unit disk with finitely many poles, and analytic at the origin. Then, there exists a hard pulse Q. 
with p rephasing time steps such that: 

. "fen* _ 

[0093] Finite rephasing time is a rather restrictive condition. Indeed Theorem 4 

indicates that for pulses with a fixed, finite rephasing time, the map Q — > Mq is essentially 
injective. By way of contrast, if one allows infinite rephasing time, then by introducing 
bound states, one can find an infinite dimensional space of pulses which produce a given 
magnetization profile. Theoretically, all but one choice of data prescribing bound states, 
with specified locations, therefore lead to pulses with infinite rephasing time. However, 
these pulses decay exponentially fast as j -> °°. Practically speaking, this means that if one 
truncates such a pulse to have a finite, but sufficiently large rephasing time, then it 
produces a very small change in the magnetization profile. The choice of data prescribing 
the bound states determines how much rephasing time is needed. Hence as a practical 
matter, there are many different pulses that produce the same magnetization profile, to 
within a specified error, and require the same number of rephasing steps. 
[0094] It is important not to confuse the magnetization profile, Mn, with the flip angle 

profile, 2arcsin \b\. In SLR, one can obtain finitely many different pulses with the same flip 
angle profile, and rephasing time, by reflecting the zeros of a across the unit circle. However, 
these pulses have different reflection coefficients, and therefore give different resulting 
magnetizations. 

CHOOSING A MAGNETIZATION PROFILE 

[0095] This section considers how to approximate the magnetization profile. At the 

outset, the ideal magnetization profile is a function defined on R. In all applications to MR, the 
profile agrees with the equilibrium state outside a finite interval. Since the invention 
approximates the magnetization profile, or data derived from it, by a periodic function, an 
interval of offset frequencies, [zq, ei], is selected that is large enough to contain all frequencies 
that arise in the application. The magnetization profile is then considered on [zo, Z\] as one 
period of a periodic function. The time step A must be chosen small enough so that z\ - Zo < A" 
l 2n. 
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[0096] Specifying a finite duration pulse. As noted above, designing a pulse, Q, of 

duration N, with p rephasing time steps, is reduced to specifying functions: 

a = 

;=o 
n ~p 

*> = Yl b J luJ ' 
/=■-/» 

with a 0 > 0, and \a\ 2 + \b\ 2 = 1. The corresponding magnetization profile is: 

f Af x 1 f 23a*b 1 

L w a J L m 2 - w* J 

[0097] In SLR, the standard method for determining polynomials a and b, which satisfy 

\a\ 2 + \b\ 2 = 1, is to first specify 1 -2\b\ 2 with \b\ < 1 as an approximation to M z . Usually, the 
approximation to b itself is a linear phase function with magnitude given by \b[ In order to 
obtain a minimum energy pulse one then sets: 

fl =exp(n+log(t-|&| 2 )). 

This resulting function a is necessarily a polynomial of degree N — 7, satisfying \a\ 2 + \b\ 2 = 
1. The polynomials a and b can then be replaced by other polynomials by reflecting some of the 
zeros across the unit circle (so that \a\ and \b\ remain unchanged). In this method, one only has 
direct control of the z-component of the resulting magnetization profile. However, direct 
control on the phase of the profile is sacrificed, in exchange for control on the duration. In 
many applications the phase errors implicitly introduced by approximating \b\ instead of b are 
close to linear within the "passband," and can therefore be reduced by rephasing. 
[0098] Specifying an infinite duration pulse with finite rephasing time. As noted 

above, in order to design a hard pulse Q with p rephasing time steps, one can specify a 

reflection coefficient re LT(S l ), with the property that rw?' 1 is meromorphic in D with 
finitely many poles, and analytic at the origin. If there are no bound states, then the Fourier 
series of r takes the form of Equation (56): 

CO 

r(ur) = £ Hj)iv>. 

The resulting magnetization is: 




One may now solve for r in terms of M y , and M z : 
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r 

V\ 

\r\ 



My + fjVf T 
" \Afy + fAf ¥ | 



\f I + A*: 



r = 

By specifying r, instead of one retains control of the entire magnetization profile. In 
fact one can use essentially the same method for approximating r as is used to obtain \b\. 
Specifying the degree of the polynomial approximating r determines the number of 
rephasing steps, but does not fix the duration of the pulse. In practice the pulses decay 
very rapidly as the index tends to — °o, and can be truncated to have essentially the same 
duration as an SLR pulse with similar design parameters. 

[0099] In the Reflection Coefficient section above, it was noted that different pulses 

can produce the same resulting magnetization. As noted below, for a given profile, the 
energy of the pulse is minimized by choosing a to be non-vanishing in D. This result, in the 
context of SLR, has been described by Pauly et al. in "Parameter Relations for the Shinnar- 
Le Roux Selective Excitation Pulse Design Algorithm," IEEE Trans, on Med. Imaging, Vol. 
10 (1991), pp. 53-65. When designing a pulse with p rephasing steps, a minimum energy 
pulse is obtained by choosing rvf' 1 to be analytic in for then a has no zeros in the unit 
disk. 

[0100] Equiripple pulse design. The SLR and 1ST approaches now will be 

compared for designing a minimum energy, selective, 90 degree pulse with a fixed 
rephasing time p in terms of the design parameters, in-slice error, out-of-slice error, and 
transition width. The ideal magnetization profile is: 





" 0 " 






I 


M < 1 




0 






= 0 = 






0 


|C| > l + T 




1 





where t > 0 is a fixed transition width. The accuracy of the two approaches may be compared using 
the following parameters: 

£i = in-slice error in M z 

$2 = out-of-slice error in \M T +iA# T |. 

For each of the two methods described below, the out-of- slice error Si is a function of the in- 
slice error S\, since x is fixed. 
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[0101] In the SLR method, one directly approximates only on the z-component of 

the ideal magnetization profile: 

0 \z\ < I 

1 id>i-HT f 



obtaining b as a linear phase function of the form 



In slice, w 1 while out of slice, b ~ 0. In practice one can find b using the Remez 
algorithm. If \M X + iM y \ < Si out of slice, then, out of slice, one needs: 



— ^ = ft*. 



Using the parameter relations from Pauly et al., one can determine the in-slice ripple S\,t in b. 
It is 

which gives the corresponding in-slice error S\. 

[0102] The second method, the 1ST method, involves designing rasa real-function of 

the form 

P -\ 

so that r ~ 1 in slice, and r ~ 0 out of slice. Again, the Remez algorithm is used. If one 
wants \M X + iM y \ < Si out of slice, then, out of slice: 

n<— % — = 

Again, using the parameter relations from Pauly et al., one can determine the in-slice ripple <5i, r 
in r as: 

I — (I -h &\ ) ^ ^ I — (I — 5*) 



which gives the corresponding in-slice error 5i 
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[0103] 



By carrying out these calculations for different values of p and x, it will be 



appreciated that the 1ST method seems to give better results (i.e. smaller values of ^1 as a 
function of Si). Figure 2 shows a comparison of in-slice versus out-of-slice errors for two 
methods of pulse design for two different choices of transition width and rephasing time. Of 
course, the first pulses technically have infinite duration though, as noted above, they can be 
truncated to have essentially the same length as the comparable SLR pulse. 
ALGORITHMS FOR GENERATING PULSES 

[0104] In this section, an algorithm for obtaining a pulse which yields a given 

magnetization profile in accordance with the invention is derived. For purposes of 
comparison, the SLR recursion is derived first. The method of the invention, referred to 
herein as the discrete inverse scattering algorithm (DIST) is then provided. The DIST is a 
particular discrete approximation of the full inverse scattering transform. As with the 
standard derivation of the SLR algorithm, these computations become far more transparent 
if the work is performed in the spin domain. The needed transformation formulae are given 
in the following section. 

[0105] The spin domain. R 3 is associated with the Lie algebra su 2 via the 

relationship: 




Then the precession operator P A corresponds to conjugation by the matrix: 




^2 j e SI/2, where w =e i; 



and the rotation corresponds to conjugation by the matrix: 



(79) 



Then: 




(80) 



More explicitly, the action of the operator PJtco is given by: 



P A R ■ -L\ A '= M,-U,l '] _i [ Af z M y -M,i] , 
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where: 

[0106] Scattering theory (part 2), £2 : z c is a hard pulse, and T_j and T+j are as 

defined above so that the scattering operator is: 

for all 7 € Z. Let us define 

The recursion in Equation (63) in terms of these operators is: 

(81) 

and the normalization is 

The scattering operator is given by 

Sa = P-j*UZ % jU-jPjA 

(82) 

for ail j eZ. 

[0107] All of these operators are now reexpressed in terms of matrices in SU2 as in the 

previous section. The scattering operator is given in terms of (a, b) by: 

The operator Sq acts on su 2 valued functions of w by conjugation by this matrix. This complex 
2x2 matrix corresponds to the real 3x3 matrix in Theorems 4 and 5. By analogy with the 
continuum case, the functions a and b are called the scattering coefficients. Similarly: 

It can be shown that, for an admissible pulse, the functions {A_j, B_j,A*+j,B*+j} extend 
analytically to the unit disk. This is the discrete analogue of the analytic extension properties 
of the normalized eigenfunctions {y/\^ V2,+ } in the continuum case. 
[0108] The recursion of Equation (81) is: 
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L 0 »" J JU' •/JL«w J I 0 W " ,/2 (83) 

where 09;, a/) are related to G$ by Equation (80). In terms of the matrix components this can 
be written as: 



where 



l *±,» J L W» » JL *±J J' 



(84) 



Notice that cycan be recovered from yby using the formula: 



V* 2\y\ 
ai = - — - a res in - 



Ik I i + IH 2 



(85) 



Inverting Equation (84) gives: 

U;HU (86) 

The scattering Equation (82) can be rewritten as: 

L m <*< J L JL *-j <; J" 

(87) 

[0109] The energy formula. Let Q. be a valid hard pulse, and let 



CO 



ft = X] ^ 

/=-co 

be the corresponding scattering coefficients. It is assumed that a has finitely many zeros in * 
{Cj} with multiplicities fdj. Then a may be written as: 

(88) 

where « is non -vanishing in ED. 
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[0110] By Equation (84): 



or 



CO 



co i CO 

(89) 



On the other hand, using Equation (88): 

m 



Since log \d\ is harmonic in the disk: 

tag™ = £ log |a I + ^- f" log \a(S 9 W0. 

/o (90) 

= £ log |ft | + ^- I** log\a{i>°y\d9 

w 27170 (91) 

Equations (89), and (91) together give: 

- J] log(l + |>v| 2 ) = -— / log latf*) \d& - 52 log |& | 

J2 log(l + IJ7I 2 ) = ^ jf o !og(l +|r(^)| 2 X/fl-252log|&| f 



r — * 

where r* ? is the reflection coefficient. This formula is the discrete analogue of the formula for 
the energy of a pulse in terms of the reflection coefficient and bound states. In the absence of 
bound states, a related, though somewhat less explicit formula appears in the aforementioned Pauly 
et al. article. 

[0111] If the hard pulse is an approximation to a smooth pulse q(t), then: 

coj ~ 2A<7(/A) 
If A is small, then Yj is also small and therefore: 

log(l + |y;i 2 )^|)';| 2 ^|^| 2 , 

so the energy of the pulse is approximated by: 



/• CJD I 6> * |2 1 C 2 ta 

|g(,)| 2 <//«- £ l-^l ~^jf >os(«+|r(c'' fl )| 2 >/fl--J]log|&|. 



(92) 
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The factors of A' 1 on the right hand side account for the scaling between the range of physical 
offset frequencies, and the interval [0, 2n]. Recall that |£k| < 1, so the inclusion of bound states 
increases the energy. 

[0112] The Shinnar-Le Roux recursion. Consider a hard pulse Q, of finite duration, so 

that a)j = Q(}) vanishes for j outside the interval [p — N, p]. The recursive part of the Shinnar-Le 
Roux (SLR) algorithm is a procedure for recovering Q from the scattering data 5q. For purposes of 
comparison, the SLR algorithm is derived using the above notation. 
[0113] Since 



r a-^k l r i i r a+. p ] 



(93) 

it can be seen from the scattering Equation (87) that: 

It follows from Equations (84) and (93) that A_j and B_j are polynomials in w for all j > p 
— N . In fact, these polynomials must be of the form: 

with A_j(0) > 0. This implies that a and b are of the form: 
with ao > 0. This proves the first part of Theorem 4. 

[0114] Using the recursion of Equation (86) together with the fact that the constant 

coefficient of B.j always vanishes: 

-a;VjA-j (0) + o,jB-j{\) = 0 ? 



which implies that: 
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One can therefore determine A.j and B.j from A_j+i and using Equation (86) along 

with this expression for Thus, one can recursively reconstruct { Yj} starting from A. iP = a 
and B iP = bvf. The pulse { a^} is then obtained using Equation (85). 

[0115] To prove Theorem 4, one still needs to show that the algorithm works if one 

starts with arbitrary functions a and b (of the correct form). The details can be found in the 
paper by Shinnar et al. entitled "Inversion of the Bloch Equation," J. Chem. Phys., Vol. 98 
(1993), pp. 6121-6128. 

[0116] The discrete inverse scattering approach. In the discrete inverse scattering 

approach of the invention, the input data is a reflection coefficient r, which is defined on the unit 
circle, along with data specifying the bound states. Half of this data is a set of points {Ci»—>Cm} 
in the unit disk. For simplicity, one may suppose that these points are distinct, though higher 
multiplicity poses no obstacle in principle, as noted below. With this simplifying assumption, 
the other half of the data is a set of m non vanishing complex numbers {ci,...,c m }. In the 
continuum case these are called norming constants. 

[0117] An alternative method for reconstructing CI from Sn will now be described that 

does not require Q to have finite duration. Given the reflection coefficient, as a Fourier series, 
this method is also recursive, though the actual recursion is quite different from that used in 
SLR. The derivation of the algorithm makes extensive usage of the projection operators, EL,, 
defined above. 

[0118] It is assumed that A_j y w' x B_j, A*+j, and B*+j are all analytic in the unit 

disk, since this is the case for admissible pulses. By the scattering Equation (87): 

r a-j "I _ r >w -n., i r * l 
L < ; J - L *+j Kj JU* ; r 

or 

: = aB+ : + hto'A* ; . 

+-i +.j (95) 

Dividing by a, rearranging, and applying n_ one gets: 



(96) 
(97) 
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h 

where r ? . The foregoing calculation is correct whether or not a vanishes in the unit disk. 
One may now derive the algorithm assuming that there are no bound states, and then give the 
necessary modifications to include this data. 

[0119] The DIST recursion with no bound states. One restricts to the case where a has 

no zeros in By analogy with the continuum theory, there are no bound states. The terms 

n -4 : i n^tt;"* 1 /?- ; 

~* •> and " rJ vanish, and so: 

These equations are the discrete analogue of the matrix Marchenko equation for the spin 
domain Bloch equation. They can be combined into the single equation: 

>*+./<*>) (9 8 ) 
By examining the coefficient of w in Equation (84): 

J (99) 

Thus jj could be obtained from the first coefficient of the solution to the Marchenko Equation 
(98). Just as in the SLR algorithm, the pulse could then be computed using Equation (85). 
One could therefore generate the pulse by solving Equation (98) directly for each j. As 
shown below, there is a recursive algorithm for determining jj. This algorithm is much more 
efficient, and stable than solving the Marchenko equation directly. 

[0120] The DIST recursion with bound states. For the case where a has zeros in D -, it 

may be assumed that there are finitely many bound states, with a vanishing on the set £ = 

{£i,...,Cm} e ^ The terms n ~° A - J and f?-*!"" 1 *-*" appearing in Equations (96) and (97) 
need to be treated. For that purpose, it is assumed that there is a sequence of functions 
{Q +J (w) : j G Z } with following properties: 

a For each j e 2 , a(w)Q+Jw) is analytic and non-vanishing in the unit disk. 

b. For each j £ the following relations hold: 

n.^.,. = -n_ e+ ,.(u,-'* +j) * (ioQ) 

a (101) 
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c For each je 2 , and k > 0 the difference Q+j+k - w k Q+j is analytic in the unit 
disk. 

It will be explained below how these functions are constructed from the data specifying the bound 
states. Since a does not vanish at w = 0, the first and third conditions show that knowledge of one 
of the functions Q+j 0 suffices to determine all the others. In practice, one such function is taken as 
the data which specifies the bound states. In the context of the present discussion, the reflection 
coefficient r(w) on the unit disk and the bound states are entirely independent of one another. In 
many applications, the reflection coefficient extends as a meromorphic function to the unit disk, and 
the data specifying the bound states is taken from the singular parts of r at these poles. 

[0121] Suppose that Q, is an admissible pulse, so that A*+j and B*+j are analytic in ° r for all 

/ With these functions, Equations (96) and (97) can be rewritten as: 

These equations can again be combined into the single Marchenko equation: 

(,+n ^ n ^S = - nw ' 

(104) 

where: 

Once again, it is possible to compute {vv,} by first solving Equation (104), and then using Equations 
(99) and (78). In practice however, the recursion derived in the next section is used. The case 

derived in the previous section, where a has no zeros in ■ can now be thought of as the special case 
where Q+j = 0 for all / 

[0122] A recursive algorithm for 1ST. A recursive algorithm is now given for the 

determination of { Equation (84) states that A+ and B+ satisfy the recursion: 

Plugging this formula into Equations (102) and (103) and using the fact that n_r 7 w = 
n_r y+ i, one gets: 

n-ry+tCui- 1 fl+j+i) 4 
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Examining the constant coefficient of the second equation: 

-VjA+J+\(0) = -»«ri<.; + |)M)-yj»(r;to«;j + |)(-l) 

As noted above, &~{h)(ri) denotes the coefficient of w" in the Fourier series of h. Solving for yj 
gives: 



(X)) - S?(r /H (11H (0) " 

(106) 



Thus, one can reconstruct A+j and 5 +J from A+j+i and B+j+\ using this expression for jj along with 
the recursion of Equation (86). To begin the recursion, it is assumed that at a sufficiently positive 
value of j, A+j is very close to 1, and B+j is very close to 0. The pulse is then obtained using 
Equation (85). 

[0123] Equation (104) is the discrete analogue of the right Marchenko equation. As in the 

continuum case, there is also a left Marchenko equation. When there are nontrivial bound states, it is 
important to use both equations in order to control the instability that arises from these ill 
conditioned linear systems. In the next section, the analogue of the left Marchenko equation is 
provided, which gives a recursion for determining { in terms of A j and B j. 
[0124] The left Marchenko equation. The same analysis is carried out beginning 

with the equations: 

Alj = *Al 9l + Vw-*B-j (1Q7) 

Blj = aBtj-b'w-U^ (10g) 

instead of Equations (94) and (95). The left Marchenko equation is then: 

(i + n + 5|n.,;> 7 £^ = --n + 5; 

A -' m (109) 

where: 

a 

I 

and Q_j is meromorphic in the disk, with the same poles and multiplicities as " The function 



Q_j is related to Q+j by: 



J J « (110) 



[0125] Using the technique from the preceding section, one obtains: 
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yf ( _ y(5^iA-. / - l )(p) 



J il _ j. i (CD *- r /-i K0>" ^ 1U ^ 

One can reconstruct { from 5 and <2_ using Equation (111) with the recursion of Equation (84). The 
recursion starts at a sufficiently negative value of j where A_j and B _j are approximated to be 1 
and 0, respectively. 
PRACTICAL PULSE SYNTHESIS 

[0126] In the previous section, recursions were derived for determining the coefficients 

{ ft], which, in turn, gives the desired hard pulse. What remains is to describe how to go from an 
ideal magnetization profile M,(z), on the real axis, and data specifying bound states, to the recursion 
described above. The data is now described for the right Marchenko equation. The calculation of the 
data for the left equation is described in the next section. 

[0127] The first step is to compute the reflection coefficient r,(z) determined by Af,(z). As 

described above, r,(z) is approximated, on a fixed finite interval, by a finite Fourier series in w = e lzA , 

(112) 

Note that, under the map z —> e lzA , the upper half plane maps to the unit disk. The time spacing, A, 
must be chosen so that the range of frequencies of interest lies in an interval of the form / = [zo - A" 
zo+A" 1 ^. The ideal reflection coefficient is assumed to vanish outside this interval. For 
minimum energy pulse synthesis, r(w) is usually found using a Remez type algorithm. For 
other applications, such as self refocused pulses, different sorts of approximations may be more 
appropriate as explained below. 

[0128] The next step is to consider how to specify the bound states. Let {zi,. . .,z m } (where 

z* = 2^) be distinct points in the upper half plane. As in the continuum case, the specification 
of the bound state data is, in principle, entirely independent of the behavior of r on the unit 
circle. If the multiplicities {dj} are all equal to one, then, by analogy with the continuum case, 
the bound states are specified by a rational function of the form: 

m 

*" ! (113) 

The location of the pole is related to Zk by = e l2 * A . The coefficients {c*} are the data 
corresponding to the norming constants {C*}. For small A, c* ~ C*A. The locations of the poles 
place a further constraint on A : the points {z/,...^ m } must lie in the strip / + iR, where I = 
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[C Vy ]. If the multiplicities are not all one, the data specifying the bound states is described by 
choosing the rational function, £>+,o, to be of the form: 

for G 5 rfjt a positive integer, and a polynomial of degree at most d k — 7, which does 
not vanish at The analogous data for the left equation is described in the next section. 
[0129] Once Q+,o is specified, then the other functions Q+j are determined by the 

relations 

Q+j = itf ; G+.o for j > 0 

Ok/io'C+.o - n ( o>(i» ; G+.o) for ; < 0. q ^ 

These functions satisfy the relations in Equations (100) and (101). Examples are given below 
that demonstrate that the inclusion of bound states dramatically alters the pulse without 
changing the reflection coefficient. In other words, the bound state data, as described herein, are 
indeed the free auxiliary parameters in the pulse synthesis problem. 

[0130] The Marchenko kernel. Using the fact that U + f = [T1J]*, it may be 

deduced that the left and right Marchenko equations only require a knowledge of Yl_rj, and 
Yl+sj, respectively. It may be assumed that rj and sj are bounded functions for all j. Since 
n_r,w* = n_r 7+ jt, whenever k > 0: 

/(«) : = »(r;)(« + j) for a + j < 0. 

This the discrete right Marchenko kernel. Similarly, the discrete left Marchenko kernel is 
defined by: 

:= S?(s;K« - /) for n - j < 0. 
[0131] Recall that r, = rw 7 " 1 - Q+j. By computing Fourier series, the right 

Marchenko kernel is given by: 

fin) = r(/i + 1) - P+jtWtT- 1 . 

where P+Jt is a polynomial of degree d k — 1. Note that the finite sum is exponentially 
decreasing as n tends to — ~ and exponentially increasing as n tends to -H>o. It is this 
exponential growth that necessitates the introduction of the left Marchenko equation when there 

are nontrivial bound states. The formula for is omitted in terms of P+jc but in the case 
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where dk = 1 it is seen that = P+* = c& a constant. The Marchenko kernel fin) may be 
nonzero for all n e 2i though in general it decays exponentially as n tends to — ©o In practice 
one truncates the kernel at a finite so that: 

CO 

22 i/c-")i 

«=A r nn+l 

is sufficiently small. As is evident in the examples, this does not present any serious 
difficulties. 

[0132] The data describing the bound states for the left Marchenko equations is 

derived from <2+,o and the scattering coefficient a. It is described by a function Q_# : 

that can be determined from Q+ t0 and a by using the fact that: 

2 

P-lP+j, = — 737- mod (u> -ft} 4 

(see Equation (110)). In the case where = 1, so P^* are constants, this becomes: 
Once Q_o is obtained, one can compute the left Marchenko kernel: 

m 

j?c«) = f(«) - 22 2 - 

Again, the formula for the polynomial is omitted except when = 1, in which case: 



p. k = = 



Thus, the role of the bound states in the structure of the Marchenko kernels exactly parallels the 
situation in the continuum theory. 

[0133] A situation will now be described where the data describing the bound states is 

derived from the reflection coefficient. Suppose that one is given an ideal reflection coefficient 
r, and would like a pulse that requires p rephasing time steps. To find a minimum energy pulse, 
r, is approximated by a function r with a Fourier series of the form: 

CO 

r(w) = 22 Hri}*f- 

;=i-/> 

This is usually done using a Remez type algorithm. One specifies transition widths, and in- 
slice and out-of-slice ripples. These parameters, along with the rephasing time, are related. In 
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general, once the rephasing time is fixed then, within certain bounds, two of the other three 
parameters can be specified. These parameter relations are described in detail by Pauly et ah 
[0134] If one wishes to include bound states, then r, is instead approximated by a 

rational function r so that rw^ 1 is analytic and nonvanishing in a neighborhood of 0. Let C, = 

{Ci>. -, Cm} be the poles in e? with multiplicities di,...,d m . It is also assumed that r does not 

vanish at for j = l,...,m. 

[0135] Define a and b as above so that r = .«* It may be shown that a and b are rational 

functions. In order for Q. to have p rephasing time steps, one sets: 

b wP~ 1 biu?' 1 

Cl+* = n-— — =n r 



a 



One can then set: 

e-., = n,__, 

because: 



b 4 b 



= n c 



a 

= n t( i- fl «) 
= n_i. 



The last step follows from the fact that a* is analytic at each (since it was assumed that r 
does not vanish at !/£**). Therefore: 

r p = 0 

s ! — n ? . 

p a * a 

is used as initial data to determine the right and left Marchenko kernels. 

[0136] As a particular case of this construction, r may be taken to be a rational 

function, analytic in a neighborhood of 0. In this case p = 1, and the corresponding pulse 
requires no rephasing, i.e., is a self refocused pulse. This is an instance where r would not be 
defined by using a Remez type algorithm. Rather one would specify the locations of the poles, 
{Ci, .., Cm}, and look for a rational function, having poles in that set, which gives an "optimal" 
approximation to r |t In general one would constrain the poles so that: 

m 

-£log|C;|<E, 
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so as to not introduce too much additional energy. This is described in detail below. 
[0137] Thus, in the case of simple bound states, the right and left kernels are related by 

a simple formula. When the multiplicities are increased, the relationship between the right and 
left kernels becomes much more complicated. However, all the information needed to derive 
the formula is contained in the above discussion. 

[0138] The explicit algorithms. Putting the pieces together, one can write down the 

explicit algorithms in terms of the left and right kernels. The reflection coefficient and bound 
state data determine the kernels: 

JBI 

fin} =?(n + \)-J2 *+Jk COtf - " 1 

and 

m 

according to the above discussion. 

[0139] One needs to select a starting point for each recursion. If fin) = 0 for rc < -no 

then Nrnax = no - 1. Otherwise the coefficients f(n) are checked and is chosen so that: 

is smaller than a fixed e > 0. Generally speaking gin) will not vanish for sufficiently 
negative n. The starting point for left recursion is therefore selected according to the latter 
criterion: N min < 0 so that: 

E l«0»>l<«. 

o=-cc 

With these choices: 

K-mvw = H 000...] 

*-+.AW = 100 00... J 

/C_.aw = II 0 00... J 

L- M < mi = 100 00... |. 

Then the following recursions are used: 

K + j(ii) = or; +■ yjL+j+t (n)J for n > 0 

L + ,j(n) = or; [-yfK+j+t (n - I) + L +J+t (« - I)] for n > I 

Z. + .;(0) = 0 

*-_.;(/!) = ft; [AT_ ri _,(n) - y?_ ,L +J for n > 0 

£._.;(») = a; [y;-|/C_.;_i(n — I) -f- L-.;-\(n - 1)] for/i > 1 
*--.;(0) = 0 
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and the formulas 

_ Y£aJt-n-j-\)K+j+H.n) 
' K+j+t (0) - E«=o /(-» -J- »)*•+. jh («) (H5) 

_ y . = S£og(-m-J-i)*-.j-iW 
(0) - E^o *(-" +7"= «)/- 

•i = ci + i^ft-^. (116) 

[0140] The point at which to switch from the right to left recursion can be determined 

by evaluating the norms of the operators defined by the Marchenko kernels. A good rule of 
thumb is to switch at an index j 0 for which: 

co jb— I 

Note that for J G Z? and n > 0: 

K^j(n) « 3^4* ; )(«) 
^ 9CA-.;)(ii) 

EXAMPLES 

Minimum energy pulses: SLR vs. 1ST. 
[0141] In these examples, both the SLR and 1ST approaches were used to design 

rather short 90° pulses with short transition regions. As a consequence, the in-slice ripple is 
fairly large. Figures 3(a) and 3(b) compare two such pulses, where the horizontal lines in the 
pulse plots, to the right of the pulse, indicate the rephasing times. Figure 3(a) illustrates SLR 
versus 1ST for a 90° pulse with 2ms rephasing time, 82 = 0.1 and the transition width is 0.2 
KHz. Figure 3(b) illustrates SLR versus 1ST for a 90° pulse with 4ms rephasing time, 62 = 
0.01 and the transition width is 0.15 KHz. They have the same rephasing times though the 
1ST pulses are overall a little longer than the SLR pulses. In both cases, the 1ST pulse produces 
cleaner transverse and longitudinal profiles. 

[0142] The effect of T2 relaxation. In this section, the effect of a nonzero spin-spin 

relaxation time on the resulting magnetization profile is studied. In computer experiments, the 

P T 

precession operator P t is replaced by the operator ' r 2 given by: 



PtS 2 M = 



[e—h 0 0 "I T 3te** 0 "I 

0 e-ri ON -%* ,z 9ta* 8 0 M. 
0 0 1 J L 0 0 U 
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Hence transverse relaxation effects are accounted for throughout the duration of the pulse. 
Figures 4(a) and 4(b) illustrate a comparison of the magnetization profiles produced by SLR 
and 1ST pulses with T 2 = 10 ms. As before, the horizontal lines in the pulse plots the rephasing 
times. Figure 4(a) illustrates T 2 = 10 ms for a 90° pulse with 2ms rephasing time, 82 = 0.1 and 
the transition width is 0.2 KHz. Figure 4(b) illustrates T 2 = 10 ms for a 90° pulse with 4ms 
rephasing time, 5 2 = 0.01 and the transition width is 0.15 KHz. Figures 4(a) and 4(b) show that 
even though the SLR pulses are overall slightly shorter, the effects of T 2 -relaxation are a little 
less pronounced on the magnetization profiles produced by the 1ST method. 
[0143] Refocusing pulse. Ideally, an inversion pulse maps the equilibrium 

magnetization [0,0,1] to the vector [0,0,-1] for offset frequencies lying in a specified band. 
For each offset frequency, such a pulse defines an orientation reversing rotation /?(£) of the 
transverse plane. The pulse is a refocusing pulse if and only if this rotation is independent of 
the offset frequency, within the band where the longitudinal magnetization is inverted. Using 
the scattering theory formulae above, it is not difficult to show that an inversion pulse is 
refocusing if and only if b(Q is real valued. This, in turn, implies that the pulse is symmetric. 
At the conclusion of such a pulse, the transformation of the transverse plane is a fixed 
reflection, e.g. 
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[0144] As noted above, it would require an infinite amount of energy to exactly invert 

the magnetization for offset frequencies belonging to a band of positive length. For this 
reason, a pulse is designed that takes [0,0,1] to [sin e,0, -cos e], for a small, positive e. It is 
not difficult to show that, if £>(£) is real, then, in-slice, this produces a family of orientation 
preserving rotations of the transverse plane of the form r € (E,)Ro, with r e (£) -Id = 0(e 2 ). 
Figure 5(a) shows a refocusing designed pulse to produce a 179°-flip of the equilibrium 
magnetization. Figures 5(b)-5(d) show the effects of this pulse on the jc, y and z components 
of the magnetization. Note the saturation that occurs outside of the slice. 
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[0145] Self-refocused pulses. In many applications it is useful to design a self- 

refocused pulse, or a pulse with no rephasing time. It is important to understand the 
relationships between pulse energy, pulse duration, and the accuracy of the approximation. For 
example, the more poles r is allowed to have, the better the approximation can be. However, 
according to Equation (92), each added pole costs energy. But perhaps the most serious issue is 
the stability of the designed pulse. As noted above, pulses designed with bound states tend to 
be more sensitive to jBi-amplitude errors. In this section, techniques to practically design 
relatively low-energy, self-refocused pulses are described. 

[0146] Let r f be the ideal reflection coefficient, as a function of the offset frequency 

z. For example, for a selective 90 degree pulse: 



n(z> = 



t if|s|<l 
0 if \z\ > 1. 



According to Theorem 4, one should approximate r, by a function r of w = e lz * which is 
meromorphic in the unit disk, and analytic at the origin. Of course, r is a periodic function of 
z, so one actually approximates r, in a finite interval containing zero. 
[0147] r is designed as a real valued rational function of w, with simple poles. It is 

awkward to do this directly, since n is a function of z. A practical alternative is to first 

approximate r, by a rational function r of the form: 



Ub-i «-«}J (li8) 

for T|i» »T|jfi in the upper half complex plane. Then: 

- V f iA< J _ ££££ . ~ l Ac 1 _ - /A *v 1 

Notice that, for a fixed z, r(e lzA ) can be made arbitrarily close to r (z) by taking A small 
enough. The contribution of the bound states to the energy of the designed pulse is then 
approximately: 

With this choice of bound state data, the right Marchenko kernel, fin) vanishes for n > 0. 

[0148] There are many possible techniques for choosing the rational function r . The 

method used to compute the examples shown in Figure 6 is to first fix m y the number of the 
poles, and then constrain the locations of the individual poles. A Newton type algorithm is then 
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used to determine the precise locations of the poles C,\ 



,. . ., so that the L 2 norm of ,r r;l can be 



made as small as possible. In general, poles close to the x-axes add relatively little energy, but 
tend to result in long pulses, whereas poles with large imaginary parts tend to yield short, high 
energy pulses. 



using different rational approximations, r, for the desired reflection coefficient. The poles of r 
in the unit disk were then used to define the bound states. As expected, a pulse was produced 
that required no rephasing time. The pulses are shown on the left of Figure 6 and the 
magnetization profiles they produce are shown on the right of Figure 6. These profiles are 
obtained without any rephasing time. 

[0150] Arbitrary profiles and bound states. A last group of examples demonstrates 

that the bound states are indeed auxiliary parameters in the pulse design problem. Figure 7 
shows four pulses, and the transverse magnetization profiles they generate. Each pulse is 
generated using a single reflection coefficient, along with data specifying bound states. The 
bound states in these examples are specified in terms of the continuum inverse scattering data, 
i.e., the location of the pole, and a norming constant, as in Equation (118). The pulses are 
evidently quite different from one another, but nonetheless, produce the same transverse 
magnetization profile. 

[0151] Figure 7 illustrates a variety of 90° pulses illustrating the effects of adding 

bound states. Figure 7(a) illustrates the minimum energy pulse with the given magnetization 
profile, M. The other three pulses were obtained by augmenting the magnetization profile data, 
with the indicated bound states. To the right of the pulse is the transverse magnetization profile 
it generates. For example, Figure 7(b) illustrates the pulse with magnetization profile, M, and a 
bound state at 0.5/, with norming constant 1.0. Figure 7(c) illustrates the pulse with 
magnetization profile, M, and bound states at 0.5/ + 1.4, and 1.0i-1.0, with norming constants 
1, and -10. Figure 7(d) illustrates the pulse with magnetization profile, M, and bound states 
and i-1, i, and i+1, with norming constants 1, 2 and 1. 

[0152] Wavelet example. In this example, a standard wavelet is used to define the de- 

magnetization. Such pulses were earlier constructed by several investigators using the SLR 
method; a discussion is provided by Panych in 'Theoretical comparison of Fourier and 
wavelet encoding in magnetic resonance imaging," IEEE Trans, on Med. Imag., Vol. 15 
(1996), pp. 141-153. If i/KS) is a wavelet, then: 



[0149] 



Figure 6 illustrates a variety of self refocused 90° pulses that were designed by 
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[0153] 



The Mexican hat wavelet is defined by the function: 




Figure 8(a) shows the minimum energy 1ST pulse producing the transverse magnetization 
defined by ^(£), which is shown in Figure 8(b). In this example, the Marchenko equations 
were solved using simple iteration. The pulse took less than a minute to compute. 
CONCLUSION 



design via the inverse scattering transform. The algorithm uses the idea of the hard pulse 
approximation, which is also employed in the SLR approach. There are two principal 
differences with the SLR approach. (1) The algorithm of the invention works, throughout the 
design process, with the full magnetization profile. In SLR only the flip angle profile is used, 
and direct control on the phase of the transverse magnetization is lost. (2) The algorithm of the 
invention uses a different recursive algorithm to recover the hard pulse from the scattering 
data. This algorithm, which is a discrete version of the Marchenko equation formalism from 
inverse scattering theory, allows for the inclusion of arbitrary bound states. The recursion in 
accordance with the algorithm of the invention is as fast as the standard SLR approach. By 
using both the left and right Marchenko equations, the algorithm of the invention controls the 
instabilities which arise in the solution of the pulse synthesis problem when bound states are 
present. 

[0155] The approach of the invention is called the discrete inverse scattering transform 

or DIST. Using the DIST, the full range of the solutions to the pulse design problem are 
recovered that are present in the continuum model. In particular, for a fixed magnetization 
profile, there is an infinite dimensional space of DIST pulses that produce the desired profile. A 
formula is provided for the energy of the pulse in terms of its magnetization profile, and the 
locations of the bound states. This formula is a direct analogue of the formula for the energy of the 
pulse in the continuum case. This formula makes explicit the dependence of the energy on the 
magnetization profile and the location of the bound states. This is very useful in the design of self 
refocused pulses with a small amount of additional energy coming from the bound states. 
[0156] In the DIST approach one does not have direct control on the duration of the pulse. 

For minimum energy pulses, however, one can control the rephasing time. In the case of minimum 
energy pulses, the effect of r 2 -relaxation is principally a function of the rephasing time. 
[0157] Thus, the present invention relates to a method for obtaining RF-pulse 

envelopes that produce a specified magnetization profile to within a predetermined error 



[0154] 



A practical algorithm is derived which implements an approach to pulse 
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based on a hard pulse approximation. The method of the invention has two facets, a 
continuum theory that describes exact solutions to idealized problems and an idealized theory 
that has a variety of different discrete manifestations which can be used to construct 
implementable algorithms. In all NMR cases, the input data is a magnetization profile. This 
is a unit 3-vector-valued function of a single real variable which, in most applications, 
vanishes outside a bounded interval. The goal of the invention is to design an RF-envelope, 
implementable on an NMR imaging device, which will a produce a good approximation to 
the desired magnetization profile. 

[0158] The continuum model provides a general solution and includes the steps of 

choosing an approximation to the magnetization profile, choosing auxiliary parameters such 
as bound states and norming constants, constructing the left and right Marchenko kernels 
from the selected data (including the auxiliary data), and solving the left and right Marchenko 
equations to give the desired RF-pulse envelope. When auxiliary parameters are present, one 
encounters significant numerical difficulties when trying to solve the Marchenko equations 
in order to obtain the pulse envelope from the magnetization profile and auxiliary data. The 
left Marchenko equation is introduced as a way to reduce the numerical instabilities. This 
leads to more stable problems, which can be solved on a digital computer. 
[0159] The analogous steps for one possible discrete implementation of this general 

solution requires an algorithm to solve the Marchenko equations. The methods of the 
invention include the use of a very efficient, stable algorithm referred to herein as the discrete 
inverse scattering transform (DIST) to solve these equations and to thereby obtain pulse 
envelopes from the magnetization profile and auxiliary parameters. The algorithm also uses a 
hard pulse approximation. 

[0160] The DIST in accordance with the invention includes the steps of 

approximating the reflection coefficient by an exponential polynomial and, if there is 
auxiliary data, encoding the auxiliary data as further exponential polynomials. Out of the 
resulting data, one constructs the left and right discrete Marchenko kernels. Using a recursive 
algorithm, derived using the inverse scattering formalism, the left and right discrete 
Marchenko equations are solved. From the solutions to the left and right discrete Marchenko 
equations, the RF-pulse envelope is constructed. 

[0161] The approach of the invention allows one to design RF-envelopes to produce 

all standard pulses used in magnetic resonance including: 
a. 9- flip for offset frequencies lying in a given band; 
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b. 0- flip for offset frequencies lying in a collections of bands; 

c. Inversion and refocusing pulses; and 

d. Minimum and maximum phase pulses. 

In each case, an important feature of the method of the invention is that the phase of the 
magnetization profile is controlled throughout the design process. The resulting signal is then 
applied to coils of an MR device for generating a desired RF field to locally modify a 
magnetic field in accordance with known techniques. 

[0162] Those skilled in the art will appreciate that there are in fact infinitely many 

solutions to the pulse design problem. Moreover, the space of solutions is described by sets of 
auxiliary parameters, and the method allows for the arbitrary specification of these 
parameters to obtain the complete set of solutions to the pulse design problem. In particular, 
the present inventors have provided a formula that relates the energy in a pulse envelope to 
the magnetization profile and the auxiliary parameters used in its design. This formula allows 
one to select the auxiliary parameters to obtain the minimum energy solution (no auxiliary 
parameters are used). 

[0163] The present invention also allows one to design pulses that cannot be obtained 

by methods that do not control the phase or allow for the specification of auxiliary 
parameters. For example, the methods of the invention may be used to design: 

a. Low energy self refocused pulses; 

b. Pairs of half pulses with an arbitrarily specified summed transverse magnetization 
profile; and 

c. Pulses where the flip angle profile is a specified wavelet, or any other function, and 
the phase is arbitrarily specified. 

Because the invention allows one to arbitrarily specify the auxiliary parameters, the methods 
of the invention may be used to obtain many self refocused pulses producing a given 
magnetization profile. Using the energy formula and constrained minimization, self- 
refocused pulses may be found that produce the desired profile, to within a specified error, 
having very little additional energy. 

[0164] Using the standard parameters used in SLR pulse design, in-slice ripple, out- 

of-slice ripple, transition width, and duration, the comparable 1ST pulses have been found to 
perform somewhat better using the techniques described herein. 

[0165] Those skilled in the art will also appreciate that numerous other modifications 

to the invention are possible within the scope of the invention. For example, the methods of 
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the invention may be applied to any two-level quantum system or subsystem described by the 
spin domain Bloch equation, such as systems in the field of spintronics and such potential 
applications as quantum computing where it is very important to control the phase of the 
excitation. Such systems may use electromagnetic waves belonging to other frequency bands 
besides RF, such as microwave or infrared, to produce the desired selective excitation. 
Accordingly, the scope of the invention is not intended to be limited to the preferred 
embodiment described above, but only by the appended claims. 
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